Methods For Subsurface Parameter Estimation In Full Wavefield Inversion And Reverse-Time Migration

ABSTRACT

Method for converting seismic data to obtain a subsurface model of, for example, bulk modulus or density. The gradient of an objective function is computed ( 103 ) using the seismic data ( 101 ) and a background subsurface medium model ( 102 ). The source and receiver illuminations are computed in the background model ( 104 ). The seismic resolution volume is computed using the velocities of the background model ( 105 ). The gradient is converted into the difference subsurface model parameters ( 106 ) using the source and receiver illumination, seismic resolution volume, and the background subsurface model. These same factors may be used to compensate seismic data migrated by reverse time migration, which can then be related to a subsurface bulk modulus model. For iterative inversion, the difference subsurface model parameters ( 106 ) are used as preconditioned gradients ( 107 ).

CROSS-REFERENCE TO RELATED APPLICATION

This application claims the benefit of U.S. Provisional Application No.61/303,148 filed Feb. 10, 2010, entitled Methods for SubsurfaceParameter Estimation in Full Wavefield Inversion and Reverse-TimeMigration, which is incorporated by reference, in its entirety, for allpurposes.

FIELD OF THE INVENTION

This invention relates generally to the field of geophysical prospectingand, more particularly, to seismic data processing. Specifically, theinvention is a method for subsurface parameter estimation in full wavefield inversion and reverse-time migration.

BACKGROUND OF THE INVENTION

Full wavefield inversion (FWI) in exploration seismic processing relieson the calculation of the gradient of an objective function with respectto the subsurface model parameters [12]. An objective function E isusually given as an L₂ norm as

$\begin{matrix}{{E = {\frac{1}{2}{\int{\int{\int{{{{p\left( {r_{g},{r_{s};t}} \right)} - {p_{b}\left( {r_{g},{r_{s};t}} \right)}}}^{2}{t}{S_{g}}{S_{s}}}}}}}},} & (1)\end{matrix}$

where p and p_(b) are the measured pressure, i.e. seismic amplitude, andthe modeled pressure in the background subsurface model at the receiverlocation r_(g) for a shot located at r_(s). In iterative inversionprocesses, the background medium is typically the medium resulting fromthe previous inversion cycle. In non-iterative inversion processes ormigrations, the background medium is typically derived usingconventional seismic processing techniques such as migration velocityanalysis. The objective function is integrated over all time t, and thesurfaces S_(g) and S_(s) that are defined by the spread of the receiversand the shots. We define K_(d)(r)=K(r)−K_(b)(r) andρ_(d)(r)=ρ(r)−ρ_(b)(r), where K(r) and ρ(r) are the true bulk modulusand density, and K_(b)(r) and ρ_(b)(r) are the bulk modulus and thedensity of the background model at the subsurface location r. We alsodefine the difference between the measured and the modeled pressure tobe p_(d)(r_(g),r_(s);t)=p (r_(g),r_(s);t)−p_(b)(r_(g),r_(s);t).

The measured pressure p, satisfies the wave equation

$\begin{matrix}{\mspace{79mu} {{{{\rho \; {\nabla{\cdot \left( {\frac{1}{\rho}{\nabla p}} \right)}}} - {\frac{\rho}{K}\overset{¨}{p}}} = {{- {q(t)}}{\delta \left( {r - r_{s}} \right)}}},\mspace{79mu} {or}}} & (2) \\{{\left( {\rho_{b} + \rho_{d}} \right)\nabla}{{{{\cdot \left( {\frac{1}{\rho_{b} + \rho_{d}}{\nabla\left( {\rho_{b} + \rho_{d}} \right)}} \right)} - {\frac{\rho_{b} + \rho_{d}}{K_{b} + K_{d}}\left( {{\overset{¨}{\rho}}_{b} + {\overset{¨}{\rho}}_{d}} \right)}} = {{- {q(t)}}{\delta \left( {r - r_{s}} \right)}}},}} & (3)\end{matrix}$

where q(t) is the source signature. By expanding the perturbation termsand keeping only the 1st order Born approximation terms, one can derivethe Born scattering equation for the pressure p_(d),

$\begin{matrix}{{{{\rho_{b}{\nabla{\cdot \left( {\frac{1}{\rho_{b}}{\nabla p_{d}}} \right)}}} - {\frac{\rho_{b}}{K_{b}}{\overset{¨}{p}}_{d}}} = {- \left\lbrack {{\frac{\rho_{b}K_{d}}{K_{b}^{2}}{\overset{¨}{p}}_{b}} - {\rho_{b}{\nabla{\cdot \left( {\frac{\rho_{d}}{\rho_{b}^{2}}{\nabla p_{b}}} \right)}}}} \right\rbrack}},} & (4)\end{matrix}$

and so p_(d) satisfies

$\begin{matrix}{{{p_{d}\left( {{r_{g}r_{s}};t} \right)} = {\int{\left\lbrack {{{\rho_{b}\left( r^{\prime} \right)}\frac{K_{d}\left( r^{\prime} \right)}{K_{b}^{2}\left( r^{\prime} \right)}{{\overset{¨}{p}}_{b}\left( {r^{\prime},{r_{s};t}} \right)}} - {{\rho_{b}\left( r^{\prime} \right)}{\nabla{\cdot \left( {\frac{\rho_{d}\left( r^{\prime} \right)}{\rho_{b}^{2}\left( r^{\prime} \right)}{\nabla{p_{b}\left( {r^{\prime},{r_{s};t}} \right)}}} \right)}}}} \right\rbrack*{g_{b}\left( {r_{g},{r^{\prime};t}} \right)}{V^{\prime}}}}},} & (5)\end{matrix}$

where V′ is the volume spanned by r′, and g_(b) is the Green's functionin the background medium.

One can derive the equations for the gradients of p_(b) using Eq. (5)and by considering the fractional change δp_(b) due to fractional changeδK_(b) and δρ_(b) over an infinitesimal volume dV,

$\begin{matrix}{{\frac{\partial p_{b}}{\partial{K_{b}(r)}} = {\frac{{\rho_{b}(r)}{dV}}{K_{b}^{2}(r)}{g_{b}\left( {r_{g},{r;t}} \right)}*{{\overset{¨}{p}}_{b}\left( {r,{r_{s};t}} \right)}}},{and}} & (6) \\{{\frac{\partial p_{b}}{\partial{\rho_{b}(r)}} = {F^{- 1}\left\{ {\frac{dV}{\rho_{b}(r)}{{\nabla{G_{b}\left( {r_{g},{r;f}} \right)}} \cdot {\nabla{P_{b}\left( {r,{r_{s};f}} \right)}}}} \right\}}},} & (7)\end{matrix}$

where P_(b)=F{x}, P_(d)=F{p_(d)}, G_(b)=F{g_(b)}, and F and F⁻¹ are theFourier transform and the inverse Fourier transform operators.

By using Eqs. 6 and 7, and using the reciprocity relationshipρ_(b)(r)G_(b)(r_(g),r)=ρ_(b)(r_(g))G_(b)(r,r_(g)),

$\begin{matrix}{\begin{matrix}{\frac{\partial E}{\partial{K_{b}(r)}} = {- {\int{\int{\int{p_{d}\frac{\partial p_{b}}{\partial{K_{b}(r)}}{t}{S_{g}}{S_{s}}}}}}}} \\{= {{- \frac{{\rho_{b}(r)}{dV}}{K_{b}^{2}(r)}}{\int{\int{\int{\left( {i\; 2\; \pi \; f} \right)^{2}{P_{b}\left( {r,{r_{s};f}} \right)}}}}}}} \\{{{G_{b}\left( {r_{g},{r;f}} \right)}{P_{d}^{*}\left( {r_{g},{r_{s};f}} \right)}{f}{S_{g}}{S_{s}}}} \\{= {\frac{{\rho_{b}(r)}{dV}}{K_{b}^{2}(r)}{\int{\int{{{\overset{.}{p}}_{b}\left( {r,{r_{s};t}} \right)}{\int{\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{g_{b}\left( {r,{r_{g};{- t}}} \right)}*}}}}}}} \\{{{{{\overset{.}{p}}_{d}\left( {r_{g},{r_{s};t}} \right)}{S_{g}}{t}{S_{s}}},}}\end{matrix}\mspace{79mu} {and}} & (8) \\\begin{matrix}{\frac{\partial E}{\partial{\rho_{b}(r)}} = {\int{\int{\int{p_{d}\frac{\partial p_{b}}{\partial{\rho_{b}(r)}}{t}{S_{g}}{S_{s}}}}}}} \\{= {{- \frac{dV}{\rho_{b}(r)}}{\int{\int{\int{{P_{d}^{*}\left( {r_{g},{r_{s};f}} \right)}{{\nabla{P_{b}\left( {r,{r_{s};f}} \right)}} \cdot}}}}}}} \\{{{\nabla{G_{b}\left( {{r_{g}r};f} \right)}}{f}{S_{g}}{S_{s}}}} \\{= {{- \frac{dV}{\rho_{b}(r)}}{\int{\int{{\nabla{p_{b}\left( {r,{r_{s};t}} \right)}} \cdot}}}}} \\{{\left\lbrack {\int{{\nabla\left( {{g_{b}\left( {r,{r_{g};{- t}}} \right)}*\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{p_{d}\left( {r_{g},{r_{s};t}} \right)}} \right)}{S_{g}}}} \right\rbrack {t}{{S_{s}}.}}}\end{matrix} & (9)\end{matrix}$

One can then use Eqs. 8 and 9 to perform full wavefield inversion in aniterative manner.

Reverse-time migration (RTM) is based on techniques similar to gradientcomputation in FWI, where the forward propagated field iscross-correlated with the time-reversed received field. By doing so, RTMovercomes limitations of ray-based migration techniques such asKirchhoff migration. In RTM, the migrated image field M at subsurfacelocation r is given as

M(r)=∫∫p _(b)(r,r _(s) ;t)∫g _(b)(r,r _(g) ;−t)*p(r _(g) ,r _(s) ;t)dS_(g) dtdS _(s),  (10)

which is very similar to the gradient equation 8 of FWI.

While Eqs. 8 and 9 provide the framework for inverting data intosubsurface models, the convergence of the inversion process often isvery slow. Also, RTM using Eq. 10 suffers weak amplitude in the deepsection due to spreading of the wavefield. Many attempts have been madeto improve the convergence of FWI or improve the amplitude of thereverse-time migration by using the Hessian of the objective function[9], i.e., a second derivative of the objective function. Computation ofthe Hessian, however, is not only prohibitively expensive incomputational resources, but it requires prohibitively large storagespace for a realistic 3-D inversion problem. Furthermore, FWI using thefull Hessian matrix may result in suboptimal inversion [2].

One may be able to perform more stable inversion by lumping non-diagonalterms of the Hessian into the diagonal terms [2]. These, however, stillrequire computation of the full Hessian matrix or at least a fewoff-diagonal terms of the Hessian matrix, which can be costlycomputations. While one may choose to use the diagonal of the Hessianonly [11], this is valid only in the high frequency asymptotic regimewith infinite aperture [1, 7].

Plessix and Mulder tried to overcome these difficulties by firstcomputing an approximate diagonal Hessian, then by scaling these byz^(α)ν_(p) ^(β), where z is the depth and ν_(p) is the compressionalwave velocity [7]. From numerical experiment, they have determined thatthe best scaling parameter is z^(0.5)ν_(p) ^(0.5). This approach,however, does not provide quantitative inversion of the subsurfacemedium parameters with correct units, since only approximate scaling hasbeen applied. Furthermore, this approach was applied to RTM where onlyvariations in compressional wave velocity is considered, and so may notbe applicable to FWI where other elastic parameters such as density andshear wave velocity vary in space.

SUMMARY OF THE INVENTION

In one embodiment, the invention is a method for determining a model ofa physical property in a subsurface region from inversion of seismicdata, acquired from a seismic survey of the subsurface region, or fromreverse time migration of seismic images from the seismic data, saidmethod comprising determining a seismic resolution volume for thephysical property and using it as a multiplicative scale factor incomputations performed on a computer to either

(a) convert a gradient of data misfit in an inversion, or

(b) compensate reverse-time migrated seismic images,

to obtain the model of the physical property or an update to an assumedmodel.

In some embodiments of the inventive method, the gradient of data misfitor the reverse time migrated seismic images are multiplied by additionalscale factors besides seismic resolution volume, wherein the additionalscale factors include a source illumination factor, a receiverillumination factor, and a background medium properties factor. Thisresults in a model of the physical property or an update to an assumedmodel having correct units.

It will be obvious to those who work in the technical field that in anypractical application of the invention, inversion or migration ofseismic data must be performed on a computer specifically programmed tocarry out that operation.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention and its advantages will be better understood byreferring to the following detailed description and the attacheddrawings in which:

FIG. 1 is a flowchart showing basic steps in one embodiment of thepresent inventive method;

FIGS. 2-5 pertain to a first example application of the presentinvention, where FIG. 2 shows the gradient of the objective functionwith respect to the bulk modulus in Pa m⁴ s, computed using Eq. 8;

FIG. 3 shows the bulk modulus update

K_(d) (r)

in Pa computed using Eq. 18 and the gradient in FIG. 2;

FIG. 4 shows the bulk modulus update

K_(d)(r)

in Pa computed using Eq. 24 and the gradient in FIG. 2;

FIG. 5 shows the gradient of the objective function with respect to thedensity in Pa² m⁷ s/kg, computed using Eq. 9;

FIGS. 6 and 7 pertain to a second example application of the presentinvention, where FIG. 6 shows the density update

ρ_(d)(r)

in kg/m³ computed using Eq. 28 and the gradient in FIG. 5; and

FIG. 7 shows the density update

ρ_(d)(r)

in kg/m³ computed using Eq. 34 and the gradient in FIG. 5.

The invention will be described in connection with example embodiments.However, to the extent that the following detailed description isspecific to a particular embodiment or a particular use of theinvention, this is intended to be illustrative only, and is not to beconstrued as limiting the scope of the invention. On the contrary, it isintended to cover all alternatives, modifications and equivalents thatmay be included within the scope of the invention, as defined by theappended claims.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

We derive inversion equations for K_(d) and ρ_(d) in the presentinvention using Eqs. 8 and 9. This is done by first taking advantage ofthe fact that p_(d) in Eqs. 8 and 9 can also be expanded using the Bornapproximation in Eq. 5. Neglecting the crosstalk components betweenK_(d) and ρ_(d), Eqs. 8 and 9 can then be approximated as

$\begin{matrix}{\frac{\partial E}{\partial{K_{b}(r)}} \approx {{- \frac{{\rho_{b}(r)}{dV}}{K_{b}^{2}(r)}}{\int{\int{\int{\int{{\frac{{\rho_{b}\left( r^{\prime} \right)}{K_{d}\left( r^{\prime} \right)}}{K_{b}^{2}\left( r^{\prime} \right)}\left\lbrack {{g_{b}\left( {r_{g},{r^{\prime};t}} \right)}*{{\overset{¨}{p}}_{b}\left( {r^{\prime},{r_{s};t}} \right)}} \right\rbrack} \times {\quad{{\left\lbrack {{g_{b}\left( {r_{g},{r;t}} \right)}*{{\overset{¨}{p}}_{b}\left( {r,{r_{s};t}} \right)}} \right\rbrack {t}{S_{g}}{S_{s}}{V^{\prime}}},\mspace{79mu} {and}}}}}}}}}} & (11) \\{\frac{\partial E}{\partial{\rho_{b}(r)}} \approx {\frac{dV}{\rho_{b}(r)}{\int{\int{\int{\int{{\frac{\rho_{d}\left( r^{\prime} \right)}{\rho_{b}\left( r^{\prime} \right)}\left\lbrack {{\nabla^{\prime}{P_{b}^{*}\left( {r^{\prime},{r_{s};f}} \right)}} \cdot {\nabla^{\prime}{G_{b}^{*}\left( {r_{g},{r^{\prime};f}} \right)}}} \right\rbrack} \times {\quad{\left\lbrack {{\nabla{P_{b}\left( {r,{r_{s};f}} \right)}} \cdot {\nabla{G_{b}\left( {r_{g},{r;f}} \right)}}} \right\rbrack {t}{S_{g}}{S_{s}}{{V^{\prime}}.}}}}}}}}}} & (12)\end{matrix}$

By changing the orders of the integral, Eq. 11 can be rewritten in thefrequency domain as

$\begin{matrix}{\frac{\partial E}{\partial{K_{b}(r)}} \approx {\frac{{\rho_{b}(r)}{dV}}{K_{b}^{2}(r)}{\int{\int{\left( {2\; \pi \; f} \right)^{4}\frac{{\rho_{b}\left( r^{\prime} \right)}{dV}\mspace{14mu} {K_{d}\left( r^{\prime} \right)}}{K_{b}^{2}\left( r^{\prime} \right)} \times {\quad{\left\{ {\int{\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{G_{b}\left( {r,{r_{g};f}} \right)}{G_{b}^{*}\left( {r_{g},{r^{\prime};f}} \right)}{S_{g}}}} \right\} \left\{ {\int{{P_{b}\left( {r,{r_{s};f}} \right)}{P_{b}^{*}\left( {r^{\prime},{r_{s};f}} \right)}{S_{s}}}} \right\} {V^{\prime}}{{f}.}}}}}}}} & (13)\end{matrix}$

The first integral term

$\begin{matrix}{\int{\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{G_{b}\left( {r,{r_{g};f}} \right)}{G_{b}^{*}\left( {r_{g},{r^{\prime};f}} \right)}{S_{g}}}} & (14)\end{matrix}$

in Eq. 13 is the approximation to the time reversal backpropagation fora field generated by an impulse source at r′, measured over the surfaceS_(g), and then backpropagated to r (See, e.g., Refs. [8, 3]). Thewavefield due to this term propagates back towards the impulse sourcelocation at r′, and behaves similar to the spatial delta functionδ(r−r′) when t=0, if the integral surface S_(g) embraces the point r′.This wavefield is correlated with the wavefield due to the second term,∫P_(b)(r,r_(s); f)P_(b)*(r′,r_(s); f)dS_(s), to form the gradient nearr=r′. The correlation of the first and the second term then decaysrapidly near r=r′. The present invention recognizes that the zone inwhich the amplitude of the correlation term is not negligible isdetermined by the seismic resolution of the survey. In the presentinvention, we make an approximation that

$\begin{matrix}{{{\int{\left( {2\pi \; f} \right)^{4}\left\{ {\int{\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{G_{b}\left( {r,{r_{g};f}} \right)}{G_{b}^{*}\left( {r_{g},{r^{\prime};f}} \right)}{S_{g}}}} \right\} \left\{ {\int{{P_{b}\left( {r,{r_{s};f}} \right)}{P_{b}^{*}\left( {r^{\prime},{r_{s};f}} \right)}{S_{s}}}} \right\} {f}}} \approx {{I_{K}(r)}{V_{K}(r)}{\delta \left( {r - r^{\prime}} \right)}}},} & (15)\end{matrix}$

where

$\begin{matrix}{{{I_{K}(r)} = {\int{\left( {2\pi \; f} \right)^{4}\left\{ {\int{\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{G_{b}\left( {r,{r_{g};f}} \right)}{G_{b}^{*}\left( {r_{g},{r;f}} \right)}{S_{g}}}} \right\} \times \left\{ {\int{{P_{b}\left( {r,{r_{s};f}} \right)}{P_{b}^{*}\left( {r,{r_{s};f}} \right)}{S_{s}}}} \right\} {f}}}},} & (16)\end{matrix}$

and V_(K) (r) is the seismic resolution at the subsurface location r.Equation 15 is equivalent to the mass lumping of the Gauss-NewtonHessian matrix by assuming that the non-diagonal components are equal tothe diagonal components when the non-diagonal components are within theseismic resolution volume of the diagonal component, and those outsidethe resolution volume zero. In other words, Eq. 15 is equivalent toimplicitly counting the number of non-diagonal components N_(i) of theGauss-Newton Hessian matrix in each i-th row that are significant inamplitude by using the seismic resolution volume of the survey, thenmultiplying the diagonal component of the i-th row by N_(i).

Seismic resolution volume may be thought of as the minimal volume at rthat a seismic imaging system can resolve under given seismic dataacquisition parameters. Two small targets that are within one seismicresolution volume of each other are usually not resolved and appear asone target in the seismic imaging system. The resolution volumes fordifferent medium parameters are different due to the difference in theradiation pattern. For example, the targets due to a bulk modulusperturbation yield a monopole radiation pattern, while those due to adensity perturbation yield a dipole radiation pattern. Seismicresolution volume V_(K)(r) for bulk modulus can be computed, forexample, using a relatively inexpensive ray approximation [6, 4].Persons who work in the technical field may know other ways to estimatethe resolution volume. For example, one may be able to empiricallyestimate the resolution volume by distributing point targets in thebackground medium, and by investigating spread of the targets in theseismic image. If the background medium contains discontinuity invelocity due to iterative nature of the inversion, the background mediummay need to be smoothed for ray tracing. One may also make a simplifyingassumption that the wavenumber coverage is uniform. The seismicresolution volume in this case is a sphere with radiusσ≈(5/18π)^(0.5)ν_(p)(r)/f_(p), where f_(p) is the peak frequency [6].One may also employ an approximation, σ≈ν_(p)(r)T/4=ν_(p)(r)/4B, where Tand B are the effective time duration and the effective bandwidth of thesource waveform, following the radar resolution equation [5].

Equation 11 can then be simplified using Eq. 15 as

$\begin{matrix}{{\frac{\partial E}{\partial{K_{b}(r)}} \approx {{- \frac{{\rho_{b}^{2}(r)}V}{K_{b}^{4}(r)}}{\langle{K_{d}(r)}\rangle}{I_{K}(r)}{V_{K}(r)}}},} & (17) \\{{and}\mspace{14mu} {so}} & \; \\{{{\langle{K_{d}(r)}\rangle} \approx {{- \frac{K_{b}^{4}(r)}{{\rho_{b}^{2}(r)}{dV}}}\frac{1}{{I_{K}(r)}{V_{K}(r)}}\frac{\partial E}{\partial{K_{b}(r)}}}},} & (18)\end{matrix}$

where

K_(d) (r)

is the spatial average of K_(d) over the seismic resolution at spatiallocation r.

Equation 16 can be further simplified if we use free-space Green'sfunction

$\begin{matrix}{{{G\left( {r_{g},{r;f}} \right)} = {\frac{1}{4\pi {{r - r_{g}}}}e^{\; k{{r - r_{g}}}}}},} & (19)\end{matrix}$

and assume that S_(g) subtends over half the solid angle. Equation 16then simplifies to

$\begin{matrix}{{{I_{K}(r)} \approx {{I_{K,s}(r)}{I_{K,g}(r)}}},{where}} & (20) \\{{{I_{K,s}(r)} = {\int{{{{\overset{¨}{p}}_{b}\left( {r,{r_{s};t}} \right)}}^{2}{t}}}},{and}} & (21) \\{{I_{K,g}(r)} = {{\int{\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{G_{b}\left( {r,{r_{g};f}} \right)}{G_{b}^{*}\left( {r_{g},{r;f}} \right)}{S_{g}}}} \approx {\frac{1}{8\pi}{\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}.}}}} & (22)\end{matrix}$

The term I_(K,s)(r) may be recognized as the source illumination in thebackground model, and I_(K,g)(r) can be understood to be the receiverillumination. One may also be able to vary the solid angle of integralat each subsurface location r following the survey geometry. Equation 11then becomes

$\begin{matrix}{\frac{\partial E}{\partial{K_{b}(r)}} \approx {{- \frac{{\rho_{b}^{2}(r)}{dV}}{K_{b}^{4}(r)}}{\langle{K_{d}(r)}\rangle}{I_{K,s}(r)}{I_{K,g}(r)}{V_{K}(r)}}} & (23) \\{and} & \; \\{{\langle{K_{d}(r)}\rangle} \approx {{- \frac{K_{b}^{4}(r)}{{\rho_{b}^{2}(r)}{dV}}}\frac{1}{{I_{K,s}(r)}{I_{K,g}(r)}{V_{K}(r)}}{\frac{\partial E}{\partial{K_{b}(r)}}.}}} & (24)\end{matrix}$

Equations 18 and 24 show that one can convert the gradient ∂E/∂K_(b)(r)into a medium parameter

K_(d)(r)

by scaling the gradient by the source and receiver illumination,resolution volume, and the background medium properties. If theinversion process is not iterative, one should be able to use Eq. 24 forparameter inversion. If the inversion process is iterative, one can use

K_(d)(r)

in Eq. 24 as a preconditioned gradient for optimization techniques suchas steepest descent, conjugate gradient (CG), or Newton CG method. It isimportant to note that Eqs. (18) and (24) yield bulk modulus with thecorrect units, i.e. are dimensionally correct, because all terms havebeen taken into account, and none have been neglected to simplify thecomputation as is the case with some published approaches. The publishedapproaches that neglect one or more of the terms source illumination,receiver illumination, background medium properties and seismicresolution volume, will not result in the correct units, and thereforewill need some sort of ad hoc fix up before they can be used foriterative or non-iterative inversion.

For density gradient, we make an assumption similar to that in Eq. 15,

∫∫[ V′P_(b)*(r′,r_(s);f)· V′G_(b)*(r_(g),r′;f)]×[ VP_(b)(r,r_(s);f)·VG_(b)(r_(g),r;f)]dS_(g)dS_(s)df≈I_(ρ)(r)V_(ρ)(r)δ(r−r′),  (25)

where

I _(ρ)(r)=∫∫∫| VP _(b)(r,r _(s) ;f)· VG _(b)(r _(g) ,r;f)|² dS _(g) dS_(s) df,  (26)

and V_(ρ)(r) is the seismic resolution for density ρ at subsurfacelocation r. The resolution volume V_(ρ)(r) differs from V_(K)(r), sincethe wavenumbers are missing when the incident and the scattered fieldare nearly orthogonal to each other. This is due to the dipole radiationpattern of the density perturbation, as discussed previously. As wasdone for V_(K)(r), one can employ ray tracing to compute the resolutionvolume V_(ρ)(r) while accounting for these missing near-orthogonalwavenumbers. Alternatively, one may be able assume V_(K)(r)≈V_(ρ)(r) byneglecting the difference in the wavenumber coverage.

The gradient equation 12 can then be rewritten as

$\begin{matrix}{{\frac{\partial E}{\partial{\rho_{b}(r)}} \approx {{- \frac{dV}{\rho_{b}^{2}(r)}}{\langle{\rho_{d}(r)}\rangle}{I_{\rho}(r)}{V_{\rho}(r)}}},} & (27) \\{{and}\mspace{14mu} {so}} & \; \\{{{\langle{\rho_{d}(r)}\rangle} \approx {{- \frac{\rho_{b}^{2}(r)}{dV}}\frac{1}{{I_{\rho}(r)}{V_{\rho}(r)}}\frac{\partial E}{\partial{\rho_{b}(r)}}}},} & (28)\end{matrix}$

where

ρ_(d)(r)

is the spatial average of ρ_(d)(r) over the seismic resolution volumeV_(ρ)(r).

We can further simplify Eq. 25 by using the vector identity(a·b)(c·d)=(a·d)(b·c)+(a×c)·(b×d) to obtain

[ V′P _(b)*(r′,r _(s) ;f)· V′G _(b)*(r _(g) ,r′;f)][ VP _(b)(r,r _(s);f)· VG _(b)(r _(g) ,r;f)]=[ V′P _(b)*(r′,r _(s) ;f)· VP _(b)(r,r _(s);f)][ V′G _(b)*(r _(g) ,r′;f)· VG _(b)(r _(g) ,r;f)]+[ V′P _(b)*(r′,r_(s) ;f)× VG _(b)(r _(g) ,r;f)]·[ V′G _(b)*(r _(g) ,r′;f)× VP _(b)(r,r_(s) ;f)].  (29)

The second term in the right-hand side of Eq. 29 is the correction termfor the dipole radiation pattern of the scattered field, and so itreaches a maximum when VP_(b) and VG_(b) are orthogonal to each other.Neglecting this correction term,

[ V′P_(b)*(r′,r_(s);f)· V′G_(b)*(r_(g),r′;f)][ VP_(b)(r,r_(s);f)·VG_(b)(r_(g),r;f)]≈[ V′P_(b)*(r′,r_(s);f)· VP_(b)(r,r_(s);f)][V′G_(b)*(r_(g),r′;f)· VG_(b)(r_(g),r;f)].  (30)

Then Eq. 12 can be approximated as

$\begin{matrix}{\frac{\partial E}{\partial{\rho_{b}(r)}} \approx {{- \frac{dV}{\rho_{b}(r)}}{\int{\int{\frac{\rho_{b}\left( r^{\prime} \right)}{\rho_{b}\left( r^{\prime} \right)}\left\{ {\int{\left\lbrack {{\nabla^{\prime}{G_{b}^{*}\left( {r_{g},{r^{\prime};f}} \right)}} \cdot {\nabla{G_{b}\left( {r_{g},{r;f}} \right)}}} \right\rbrack {S_{g}}}} \right\} \times \left\{ {\int{\left\lbrack {{\nabla^{\prime}{P_{b}^{*}\left( {r^{\prime},{r_{s};f}} \right)}} \cdot {\nabla{P_{b}\left( {r,{r_{s};f}} \right)}}} \right\rbrack {S_{s}}}} \right\} {f}{V^{\prime}}}}}}} & (31)\end{matrix}$

We can approximate the integral over S_(g) in Eq. 26 by using free-spaceGreen's function,

$\begin{matrix}{{{\int{{{\nabla{G_{b}^{*}\left( {r_{g},{r;f}} \right)}} \cdot {\nabla{G_{b}\left( {r_{g},{r;f}} \right)}}}{S_{g}}}} \approx {\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{\int{{{\nabla{G_{b}^{*}\left( {r_{g},{r;f}} \right)}} \cdot {\nabla{G_{b}\left( {r,{r_{g};f}} \right)}}}{S_{g}}}}} \approx {\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}\frac{\left( {2\pi \; f} \right)^{2}}{8\pi \; {v_{p}^{2}(r)}}}},} & (32)\end{matrix}$

The integral was performed over half the solid angle under theassumption that V(ρ_(b)(r_(g))/ρ_(b)(r))≈0, and ρ_(b) (r_(g)) isconstant along S_(g).

The gradient equation above can then be rewritten as

$\begin{matrix}{\frac{\partial E}{\partial{\rho_{b}(r)}} \approx {{- \frac{dV}{\rho_{b}^{2}(r)}}{\langle{\rho_{d}(r)}\rangle}{I_{\rho,s}(r)}{I_{\rho,g}(r)}{V_{\rho}(r)}}} & (33) \\{{and}\mspace{14mu} {so}} & \; \\{{{\langle{\rho_{d}(r)}\rangle} \approx {{- \frac{\rho_{b}^{2}(r)}{dV}}\frac{1}{{I_{\rho,s}(r)}{I_{\rho,g}(r)}{V_{\rho}(r)}}\frac{\partial E}{\partial{\rho_{b}(r)}}}},} & (34) \\{where} & \; \\{{{I_{\rho,s}(r)} = {\int{{{\nabla{\overset{.}{p}\left( {r,{r_{s};t}} \right)}}}^{2}{t}{S_{s}}}}},} & (35) \\{and} & \; \\{{I_{\rho,s}(r)} = {\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{\frac{1}{8\pi \; {v_{p}^{2}(r)}}.}}} & (36)\end{matrix}$

As was the case with

K_(d)(r)

, Eqs. 28 or 34 can be used as an inversion formula for non-iterativeinversion, or as a preconditioned gradient equation for iterativeinversion. It is important to note that these equations yield densitywith the correct units, i.e. are dimensionally correct, because allterms have been taken into account, and none have been neglected tosimplify the computation as is the case with some published approaches.The same is true for Eqs. 18 and 24 for bulk modulus. The publishedapproaches that neglect one or more of the terms source illumination,receiver illumination, background medium properties and seismicresolution volume, will not result in the correct units, and thereforewill need some sort of ad hoc fix up before they can be used foriterative or non-iterative inversion.

Since the first iteration of FWI is similar to RTM, the method providedhere can be applied to analyze the amplitude term in RTM with littlemodification. Seismic migration including RTM is typically used to imagethe structure of the subsurface, and so amplitude information in themigrated image is often discarded. We show that the RTM amplitude, whenproperly scaled using the method provided here, represents thedifference between the true compressional wave velocity of thesubsurface and the velocity of the background model.

We note that the RTM equation 10 is missing the double derivative of theincident field in Eq. 8. This double derivative represents that highfrequency components scatter more efficiently than the low frequencycomponents in the classical Rayleigh scattering regime (Refs. [10, 13]).Therefore, we can consider Eq. 10 as the gradient computation operationin Eq. 8, partially neglecting the frequency dependence of the scatteredfield,

$\begin{matrix}{{\frac{\partial E}{\partial{K_{b}(r)}} \approx {{- \frac{\rho_{b}{{dV}\left( {\; 2\pi \; f_{c}} \right)}^{2}}{K_{b}^{2}(r)}}{\int{\int{\int{{P_{b}\left( {r,{r_{s};f}} \right)}{G_{b}\left( {r,{r_{g};f}} \right)}{P^{*}\left( {r_{g},{r_{s};f}} \right)}{f}{S_{g}}{S_{s}}}}}}} \approx {{- \frac{\rho_{b}{{dV}\left( {\; 2\pi \; f_{c}} \right)}^{2}}{K_{b}^{2}(r)}}{M(r)}}},} & (37)\end{matrix}$

where f_(c) is the center frequency of the source waveform. Frequencydependence is partially neglected because, while the frequencydependence in the forward field p_(b) has been neglected, the frequencydependence implicit in the received field p_(s) cannot be neglected.Spatial variation of density is usually not considered in RTM, and soρ_(b) is assumed to be constant in Eq. 37.

One can now apply the same approximation used to derive Eqs. 17 and 23to Eq. 10,

$\begin{matrix}{{\frac{\partial E}{\partial{K_{b}(r)}} \approx {{- \frac{\rho_{b}^{2}{I_{K,g}(r)}{V_{K}(r)}{dV}{\langle{K_{d}(r)}\rangle}\left( {{2\pi}\; f_{c}} \right)^{2}}{K_{b}^{4}(r)}}{\int{\int{\left( {\; 2\pi \; f} \right)^{2}{P_{b}\left( {r,{r_{s};f}} \right)}{P_{b}^{*}\left( {r,{r_{s};f}} \right)}{f}{S_{s}}}}}}},} & (38)\end{matrix}$

which, together with Eq. 37, yields

$\begin{matrix}{{\langle{K_{d}(r)}\rangle} \approx {\frac{K_{b}^{4}(r)}{\rho_{b}{I_{K,g}(r)}{V_{K}(r)}}{\frac{M(r)}{\int{\int{{{{\overset{.}{p}}_{b}\left( {r,{r_{s};t}} \right)}}^{2}{t}{S_{s}}}}}.}}} & (39)\end{matrix}$

Equation 39 enables quantitative analysis of the amplitude in thereverse-time migrated image. More specifically, it enables inversion ofthe amplitude into the difference bulk modulus of the subsurface.

FIG. 1 is a flowchart showing basic steps in one embodiment of thepresent inventive method. In step 103, the gradient of an objectivefunction is computed using an input seismic record (101) and informationabout the background subsurface medium (102). At step 104, the sourceand receiver illumination in the background model is computed. At step105, the seismic resolution volume is computed using the velocities ofthe background model. At step 106, the gradient from step 103 isconverted into the difference subsurface model parameters using thesource and receiver illuminations from step 104, seismic resolutionvolume from step 105, and the background subsurface model (102). If aniterative inversion process is to be performed, at step 107 thedifference subsurface model parameters from step 106 are used aspreconditioned gradients for the iterative inversion process.

EXAMPLES

We consider the case of a 30 m×30 m×30 m “perfect” Born scatterer in ahomogeneous medium with K_(b)=9 MPa and ρ_(b)=1000 kg/m³. The target iscentered at (x,y,z)=(0,0,250 m), where x and y are the two horizontalcoordinates and z is the depth. The target may be seen in FIGS. 2-7 as a3×3 array of small squares located in the center of each diagram. Weassume co-located sources and receivers in the −500 m≦x≦500 m and −500m≦y≦500 m interval with the 10-m spacing in both x and y directions. Weassume that the source wavelet has a uniform amplitude of 1 Pa/Hz at 1 min the 1 to 51 Hz frequency band.

In the first example, we assume that the target has the bulk modulusperturbation given as K_(d)=900 kPa. FIG. 2 shows the gradient∂E/∂K_(b)(r) along the y=0 plane using Eq. 8. The scattered field p_(d)in Eq. 8 has been computed using Eq. 5. The gradient in FIG. 2 has unitsof Pa m⁴ s, and so cannot be directly related to K_(d). As described inthe “Background” section above, this is a difficulty encountered in someof the published attempts to compute a model update from the gradient ofthe objective function.

FIG. 3 shows

K_(d)(r)

using Eq. 18 of the present invention. It has been assumed that theseismic resolution volume V_(K)(r) is a sphere with radiusσ=ν_(p)(r)/4B=15 m. One can see that FIG. 3 is the smeared image of thetarget, since

K_(d)(r)

is the averaged property over the seismic resolution volume. FIG. 4 is

K_(d) (r)

when the less rigorous Eq. 24 is employed. One can see that

K_(d)(r)

in FIGS. 3 and 4 are in good agreement with each other. The value of

K_(d)(r)

at the center of the target in FIG. 3 is 752 kPa, and that in FIG. 4 is735 kPa, both of which are within 20% of the true value of 900 kPa.

The second example is the case where the target has a densityperturbation of ρ_(d)=100 kg/m³. FIG. 5 shows the gradient ∂E/∂ρ_(b)(r)along the y=0 plane using Eq. 9. The scattered field in Eq. 9 has beencomputed using Eq. 5. The gradient in FIG. 5 has units of Pa² m⁷ s/kg.As with the first example, the different units prevent the gradient frombeing directly related to a density update, again illustrating theproblem encountered in published methods.

FIG. 6 shows

ρ_(d)(r)

using Eq. 28. The seismic resolution volume V_(ρ)(r) has been assumed tobe identical to V_(K)(r). FIG. 7 is

ρ_(d)(r)

when Eq. 34 is employed. Estimation of

ρ_(d)(r)

using Eq. 34 results in less accurate inversion than that using Eq. 28due to the neglected dipole illumination term in Eq. 31.

The foregoing patent application is directed to particular embodimentsof the present invention for the purpose of illustrating it. It will beapparent, however, to one skilled in the art, that many modificationsand variations to the embodiments described herein are possible. Allsuch modifications and variations are intended to be within the scope ofthe present invention, as defined in the appended claims. Personsskilled in the art will readily recognize that in practical applicationsof the invention, at least some of the steps in the present inventivemethod are performed on or with the aid of a computer, i.e. theinvention is computer implemented.

REFERENCES

-   [1] G. Beylkin, “Imaging of discontinuities in the inverse    scattering problem by inversion of a causal generalized Radon    transform,” J. Math. Phys. 26, 99-108 (1985).-   [2] G. Chavent and R.-E. Plessix, “An optimal true-amplitude    least-squares prestack depth-migration operator,” Geophysics 64(2),    508-515 (1999).-   [3] D. R. Jackson and D. R. Dowling, “Phase conjugation in    underwater acoustics,” J. Acoust. Soc. Am. 89(1), 171-181 (1991).-   [4] I. Lecomte, “Resolution and illumination analyses in PSDM: A    ray-based approach,” The Leading Edge, pages 650-663 (May 2008).-   [5] N. Levanon, Radar Principles, chapter 1, pages 1-18, John Wiley    & Sons, New York (1988).-   [6] M. A. Meier and P. J. Lee, “Converted wave resolution,”    Geophysics 74(2), Q1-Q16 (2009).-   [7] R. E. Plessix and W. A. Mulder, “Frequency-domain    finite-difference amplitude-preserving migration,” Geophys. J. Int.    157, 975-987 (2004).-   [8] R. P. Porter, “Generalized holography with application to    inverse scattering and inverse source problems,” Progress in Optics    XXVII, E. Wolf, editor, pages 317-397, Elsevier (1989).-   [9] R. G. Pratt, C. Shin, and G. J. Hicks, “Gauss-Newton and full    Newton methods in frequency-space seismic waveform inversion,”    Geophys. J. Int. 133, 341-362 (1998).-   [10] J. W. S. Rayleigh, “On the transmission of light through an    atmosphere containing small particles in suspension, and on the    origin of the blue of the sky,” Phil. Mag. 47, 375-384 (1899).-   [11] C. Shin, S. Jang, and D.-J. Min, “Waveform inversion using a    logarithmic wavefield,” Geophysics 49, 592-606 (2001).-   [12] A. Tarantola, “Inversion of seismic reflection data in the    acoustic approximation,” Geophysics 49, 1259-1266 (1984).-   [13] R. J. Urick, Principles of Underwater Sound, chapter 9, pages    291-327, McGraw-Hill, New York, 3rd edition (1983).

1. A method for determining a model of a physical property in asubsurface region from inversion of seismic data, acquired from aseismic survey of the subsurface region, or from reverse time migrationof seismic images from the seismic data, said method comprisingdetermining a seismic resolution volume for the physical property andusing it as a multiplicative scale factor in computations performed on acomputer to either convert a gradient of data misfit in an inversion, orcompensate reverse-time migrated seismic images, to obtain the model ofthe physical property or an update to an assumed model.
 2. The method ofclaim 1, further comprising multiplying the gradient of data misfit orthe reverse time migrated seismic images by additional scale factorsincluding a source illumination factor, a receiver illumination factor,and a background medium properties factor to then obtain the model ofthe physical property or an update to an assumed model.
 3. The method ofclaim 1, wherein the seismic resolution volume is determined by raytracing using velocity of a background medium model and using an assumedfunction of frequency for a seismic wavelet.
 4. The method of claim 2,wherein the model is determined from inversion of seismic data, saidmethod further comprising: assuming an initial model of the subsurfaceregion specifying a model parameter at discrete cell locations in thesubsurface region; forming a mathematical objective function to measuremisfit between measured seismic data and model-calculated seismic data;selecting a mathematical relationship that gives an adjustment, i.e.update, to the initial model that would reduce the misfit, saidmathematical relationship relating said adjustment to a scaled gradientof the objective function, said gradient being with respect to saidmodel parameter, the scaling comprising the four scale factors, i.e. theseismic resolution volume factor, the source illumination factor, thereceiver illumination factor, and the background medium propertiesfactor, all of which appear in the mathematical relationship asmultiplicative factors that scale the gradient of the objective functionto yield the adjustment to the model parameter; and using a computer tocompute the adjustment from the mathematical relationship, and thenupdating the initial model with the computed adjustment.
 5. The methodof claim 4, wherein the physical property, i.e. the model parameter, iseither bulk modulus or density, or a combination of bulk modulus anddensity.
 6. The method of claim 4, wherein the mathematical relationshipdepends upon the physical property.
 7. The method of claim 4, whereinthe background medium properties factor comprises bulk modulus raised tothe fourth power divided by density squared when the physical propertyis bulk modulus, and comprises density squared when the physicalproperty is density.
 8. The method of claim 4, wherein the receiverillumination factor when the physical property is bulk modulus isapproximated by (1/8π)(ρ_(b)(r_(g))/ρ_(b)(r)), where ρ_(b)(r) isbackground density at location r and r_(g) is the receiver's location;and the receiver illumination factor when the physical property isdensity, I_(ρ,g)(r), is approximated by:${I_{\rho,g}(r)} = {\frac{\rho_{b}\left( r_{g} \right)}{\rho_{b}(r)}{\frac{1}{8\pi \; {v_{p}^{2}(r)}}.}}$where ν_(p)(r) is velocity at location r.
 9. The method of claim 4,further comprising repeating the method for at least one iteration,where in the initial model is replaced by the updated model from theprevious iteration.
 10. The method of claim 9, wherein the objectivefunction's functional form and the mathematical relationship'sfunctional form are unchanged from one iteration to a next iteration.11. The method of claim 4, wherein the adjustment to the initial modelis computed by minimizing the objective function using the objectivefunction's Hessian resulting in an Hessian matrix, wherein off-diagonalelements of the Hessian matrix are ignored when outside of a seismicresolution volume.
 12. The method of claim 11, wherein off-diagonalelements of the Hessian matrix within a seismic resolution volume areassumed to be equal to a corresponding diagonal element, resulting incomputing only diagonal elements.
 13. The method of claim 2, wherein thereceiver illumination factor when the physical property is bulk modulusis approximated by an integral of a Green's function in free space overa surface defined by the seismic survey's receiver spread.
 14. Themethod of claim 4, wherein the receiver illumination factor when thephysical property is density is approximated by an integral of agradient of a Green's function in free space over a surface defined bythe seismic survey's receiver spread.
 15. The method of claim 1, whereinthe seismic resolution volume is approximated as a sphere based on anassumption of uniform wavenumber coverage.
 16. The method of claim 4,further comprising using the updated model to precondition the gradientin an iterative optimization technique.
 17. The method of claim 2,wherein the model is determined from reverse time migration of seismicimages from the seismic data.
 18. The method of claim 17, wherein thephysical property is bulk modulus.
 19. The method of claim 17, whereinthe background medium properties factor comprises bulk modulus squareddivided by density.
 20. The method of claim 17, wherein the receiverillumination factor is approximated by (1/8π)(ρ_(b)(r_(g))/ρ_(b)(r)),where ρ_(b)(r) is background density at location r and r_(g) is thereceiver's location.
 21. The method of claim 17, wherein the migratedimages seismic amplitudes are converted into difference bulk modulus ordifference compressional wave velocity.